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ABSTRACT 

<D 

\ For the first time, we report a unified microscopic-macroscopic Monte Carlo sim- 

ulation of gas-grain chemistry in cold interstellar clouds in which both the gas-phase 
and the grain surface chemistry are simulated by a stochastic technique. The surface 
chemistry is simulated with a microscopic Monte Carlo method in which the chemistry 
occurs on an initially flat surface. The surface chemical network consists of 29 reactions 
initiated by the accreting species H, O, C, and CO. Four different models are run with 
diverse but homogeneous physical conditions including temperature, gas density, and 
diffusion-barrier-to-desorption energy ratio. As time increases, icy interstellar mantles 
begin to grow. Our approach allows us to determine the morphology of the ice, layer 
by layer, as a function of time, and to ascertain the environment or environments for 
individual molecules. Our calculated abundances can be compared with observations of 
ices and gas-phase species, as well as the results of other models. 

Subject headings: ISM: clouds, ISM: molecules, ISM: molecular processes 



1. Introduction 

Interstellar molecular clouds are uniquely large chemical systems. For instance, giant molecular 
clouds are typically of dimension 20-100 pc while molecular cores have a size ~ 0.1 pc. Yet, in 
these large entities, species such as H2 and ices such as water are mainly formed on surfaces of tiny 
dust particles with radii ranging in size from less than O.Ol^t to somewhat more than 0.1//. Species 
are exchanged among numerous independent dust grains and a large gas-phase body by accretion 
and desorption processes. 

The extent of a chemical reaction system in terms of density and physical size can determine 
the mathematical approach by which we calculate the temporal evolution of chemical abundances. 
Chemical reaction events are essentially discrete processes, so fluctuations exist for the evolution of 
chemical species in any system. Because molecular clouds contain large amounts of gaseous material, 
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however, the fluctuations due to randomness of the chemical reaction events in the gas phase cancel 
out, so that we can adopt the rate -equation (RE) approach to calculate the evolution of species in 
the gas phase (jGarrod et al.ll2008l ). One might think as well that because there are so many dust 
particles in clouds, the fluctuations of chemical species on each dust grain can also cancel out so 
rate equations are still applicable. It turns out that this hypothesis is true only when the averag e 
number of individual reactive species per dust grain is not well below 1 (jTielens Charnlevlll997l ). 
When the average number of reactive species is significantly below unity, significant mistakes can be 



made with the RE approach (Tielens & Charnlcy 



1997; 



20031 : iLipshtat & Bihamlhooilstantcheva k Herbstll2004l ) 



Caserii et al.lll998l ; iHerbst k Shematovich 



Various approaches have been suggested and tested to solve this problem of the surface chem- 
istry on dust particles. The modified rate equation approach is t he most computationally efficient 



one and can be imp lemented in a large chemical reaction network (jCaselli et al.lll998l : iGarrod 1120081 : 
Garrod et al.ll2009i ). However, because of its empirical nature, the correctness of this approach can- 



not be guaranteed in all instances, and should be tested by more rigorous approaches, which are 
stochastic in nature. Solution of the chemical master equation, which describes the evolution of 
probability of so-called occupancy states, is much more rigorous. There are two ways to solve the 
chemical master equation, and both are computationally tedious. The so-called "macroscopic Monte 
Carlo approach" was the first to be applied in a strochemical simulati ons. Initially, it was used to 
simulate gas-phase or surface kinetics separately (jCharnleyl ll998. 



2001 



, but eventu ally was used to 



simu l ate the chemistry occurrin g with a full gas-grain chemical reaction network (jVasyunin et al 



2009 



Charnley k Rodgerd 12009). The advantage of a Monte Carlo simulation is that it actually is 



a numerical experiment that realizes one possible trajectory for the evolution of occupancy states 
of different species; thus, one does not have to c onsider all possible occupancy states. The chemical 



master equation can also be directly integrated ( Biham et al. 2001 ; Green et alJboOll ). The advan- 



tage of direct integration is that it can be easily coupled to RE simulation of gas-phase kinetics. 
Because of the large number of possible occupancy states, how ever, approximations have t o be made 
and so far this method has only been used in small systems (jStantcheva k Herbstl 12004 ). There is 
another approximate app roach, the moment equation method, which, at low order, is much more 
computationally efficient (jLipshtat k Bihamll2003l : iBarzel k Biham 1120071 ) . but is also less rigorous. 
Recently, it has been successfully applied in a hybrid approach in which it is turned on and off 
dependin g upon the numbe r of species per grain to simulate a large gas-grain chemical reaction 



network (IDu k Parisdl201ll ) . The lowest order moment equation approach bears a strong similarity 



to the most recent modified rate equation approach (A. I. Vasyunin 2012, private communication). 

Another problem that occurs in the simulation of gas-grain chemistry derives from the com- 
plexity of surfaces. Surfaces can be rough; on such surfaces, diffu sion and desorption rat es for 
species can vary depending on the local position of each particle (jCuppen k Herbstl 120051 ) . Al- 
though the diffusive, or Langmuir-Hinshelwood, mechanism for surface chemistry is the dominant 
and best understood process, other mechanisms exist, including the Eley-Rideal process, in which 
a gas-phase species lands directly on an adsorbate, and the hot-atom mechanism, in which a gas- 
phase species lands obliquely on a surface and travels in a non-thermal manner for a considerable 
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distance. So far, the only method that can rigorously treat problems with complex surfaces is 
the continuous-time random walk (CTRW) Monte Carlo simulation approach (jChang et al.ll2005l ). 
also known variously as a microscopic Monte Carlo or kinetic Monte Carlo approach. With the 
microscopic Monte Carlo approach, one can keep track of the position of each particle on a surface, 
which is represented as a planar lattice, so that processes such as diffusion and desorption can be 
tracked based on the local environment of each adsorbate. The microscopic Monte Carlo simu- 



lation was initially used to si mulate molecular hydrogen formati on on dust grains (jChang et al 



20051 ; ICuppen k, Herbstl 120051 ) . and is still used for this purpose (jlqbal et al.ll2012l ). Inclusion of 
rough or amorphous surfaces extends the t emperature range o f efficient Ho form ation at low tem- 
peratures, in agreement with experiment (jChang et al.ll2005l ; IVidali et al.l 12009). The technique 
has also been applied to simulate surface chemistry occurring with larger surface reaction networks 
(jCuppen et al.ll2009l : ICazaux et al.ll2010l ). Based on previous simulations, it would appear that the 
complexity of surfaces cannot necessarily be ignored. 

Up to now, most Monte Carlo simulations of interstellar chemistry have contained the assump- 
tion that the abundances of gas-phase species stay constant instead of simulating a full gas-grain 
reaction network. It is difficult to couple a microscopic or macroscopic Monte Carlo simulation 
of surface reactions with an RE simulation of reactions in the gas phase because the RE simu- 
lation describes an infinitely large system without fluctuations while the Monte Carlo simulation 
describes a finite system with fluctuations. Mathematically, the Monte Carlo simulation gives the 
abundances of species as discrete numbers while the abundances of species by an RE simulation 
are continuous real numbers. Neverthele ss, we developed a n iterative method that combines these 
two approaches for the microscopic case (jChang et al.ll2007l ). The approach is successful only when 
the ratio of the accretion rate to the desorption rate of species remains constant or the desorption 
rate of species from grain surfaces is negligible. An alter native is to use a Mon te Carlo approach 
for both the surface and gas-phase chemistries. Recently IVasvunin et al.l (|2009l ) succeeded in such 
an approach using a macroscopic Monte Carlo method; they isolated a finite volume of gas phase 
species surrounding a dust particle and did a full gas-grain Monte Carlo simulation. (See also 
Charnlev h Rodgersl (|2009l ).) The agreement with the RE approach depends upon the physical 
parameters assumed; for conditions in which a stochastic treatment is needed, dramatic differences 
in abundances can occur. Closer agreement with the results of a modified rate method for the 
surface chemistry co upled with an RE treatment for the gas chemistry was subsequently obtained 
(|Oarrod et al.ll2009h . 



In this paper, we report the first successful unified microscopic-macroscopic Monte Carlo sim- 
ulation with a gas-grain reaction network. In this approach, gas-phase kinetics are treated macro- 
scopically while grain-surface kinetics on a rough surface are treated by a microscopic simulation. 
The gas- grain model that we used is introduced in Section 2, while our Monte Carlo technique is 
introduced in Section 3. Section 4 contains the results of our simulation, and a comparison with 
results from other techniques as well as observations is reported in Section 5. Section 6 contains a 
discussion along with our conclusions. 
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2. Gas-grain Model 



Gas phase species and reactions are taken from the benchmark study of lSemenov et al.l (|2010l ). 
which contains 4406 reactions and 459 species. We choose physical conditions to pertain to cold 
dense cloud cores, with a standard density of hydrogen nuclei nn of 2 x 10 4 cm~ 3 , and temperatures 
between 10 K and 16.5 K because some parameters are available only within this temperature 
range, as explained later. Each model is run at a fixed temperature; we do not consider warm-up 
and possible ice evaporation due to proximity to protostellar sources of radiation. We maintain 
a standard dust-to gas number density ratio of 10~ 12 , with dust grains of radius O.lfi. We vary 
nn to calculate the dependence of the abundances and environments of surface species on this 
parameter. The visual extinction, A v , is fixed at 10 mag, while the cosmic ray ioniz ation rate Cho is 
1.3 x 1 0~ 17 s _1 . Initial gas-phase abundances are given in Table [TJ these, taken from lSemenov et al 



(120101) . refer to low-metal ab undances except for values based on new measurements of C, N, and 
O dWakelam fc Herbstl 120081 ) . 

We start from a flat surface and choose the density of sites to be s = 1.5 x 10 15 cm -2 , as 
measured for olivine ( Katz et al.lll999l ). Thus there are 1.88 x 10 6 sites per dust grain. We assume 
that the flat surface consists of unit surface cells such that each site has four nearest neighbors. As 
the ice builds up, we assume a simple cubic structure in which each species can have six nearest 
neighbors, including one above and one below the plan e in which the spec ies resides. The detailed 
surface model has been explained in our earlier paper (jChang et al.1 120071 ). The desorption energy 
and diffusion barrier are the main parameters that determine the mobility of species on grain 
surfaces. On a rough surface, these parameters can vary with position, so that mobility depends 
on the local environment. The values of the desorption energy and diffusion barrier come from two 
contributions. The vertical interaction is that between the adsorbate and the species lying below 
it. We use E^ and EP to represent the vertical desorption energy and vertical diffusion barrier, 
respectively. The other contribution derives from lateral bonds, which are from interactions between 
the species and their horizontal neighbor molecules. Each molecule on a site can form up to to 4 
lateral bonds with its neighbor molecules. The actual number of bonds formed is dependent on the 
location of the molecule. The strength of lateral bonds is not known very well. In order to explain 
the high efficiency of H2 formation over a wide temperature range, the strengt h of lateral bonds 
cann ot be too small in order to trap enough H atoms on the grain surface (jCuppen &: Herbst 



2005). However, if the lateral bond strength is very large, species on a grain surface will lose 



much of their mobility, and so-called "islands" of adsorbates will tend to form, making the surface 
effectively "rough." We choose E^, the lateral bond strength for any species i, t o be 10% of its 
vertical desorption energy, E D , a minimal value according to ICuppen Herbstl (|2005l ) , but one 
that still broadens t he temperature range for efficient H2 formation. The local desorption energy 
is Ed = E^ + n£"iL ( Cuppen Herbstl 12005). where n is the number of lateral bonds formed by 
the s pecies and neighb o r mol ecules. Similarly, E\> = E® + nE^. Table [2] gives E D values taken 
from iGarrod Herbstl (|2006l ) which we use for each surface species considered here. We choose 



p = E®/E D to be 0.77 or 0.5 in our simulations, which corresponds to a slow diffusion rate and a 
fast rate, respectively, assuming the better known desorption energy to be fixed. The high value 
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derives from the experiment on H2 formation on olivi ne by iKatz et al.l ()1999l ) , whereas the low 
value is more suitable for ices (jGarrod HerbstJ 120061 ) . The rate that a species hops from one 
site to a nearest neighbor site is given by the expression b\ = z^exp(— Eb/T), where Eb is the local 
diffusion barrier, T is the grain temperature and v is the trial frequency, which is normally chosen 
to be 10 12 s _1 for physisorbed species. Similarly, the rate of thermal desorption is expressed by the 
relation 62 = ^exp(— E-q/T). We assume there is no tunneling under the diffusi on barrier becaus e 
none was found for analogous systems (e.g. olivine) studied in the laboratory (|Katz et al.l 11999) . 
We also assume that the products of surface reactions other than H2 do not leave the surface with 
energy taken from the exot hermicity of reaction, a lthough it is possible that they do s o according 



to a s tatistical mechanism (jGarrod &: Herbstl 120061 ) . Based partially on the analysis of IKatz et al 



(|1999l ) and partially on the need to minimize the simulation time, however, we allow all newly 
formed H2 to leave the surface. 

Photodesorption is included in our simulation. The flux of FUV photons is given by 

F FUV = G F e^ Av + G' F , (1) 

where Fq = 10 8 photons cm~ 2 s™ 1 represents the standard interstellar radiation field, while Go and 
Gq are scal ing factors. We ado pt the value for the exponential factor 7 = 2 from the experimental 
analysis of lOberg et al.l (|2007l ). The first term in the equation is for ext ernal photons whil e the 
second term is for cosmic ray induced photons. Based on a calculation by IShen et al. I (|2004l ). the 
scaling factor G for cosmic ray induced photons is 4 orders of magnitude smaller than the external 
radiation scaling factor Go = 1. 

Table [3] shows the chemical reaction network that is used for grain-surface chemistry in our 
model. The reactions are of two types: those needed to produce major species such as water, 
methanol, and methane and used previously, and those needed to eliminate high abundances of 
radicals, which can otherwise get trapped in the ice, especially at low temperatures. The destruction 
of radicals is assumed to lead to non-radical products, as found in gas-phase databases, rather than 
to other radicals. For example, OH and CH2 do not add to form the radical CH2OH but produce 
H2CO and H. Although the addition reactions to form more complex radicals are doubtless present, 
their exclusion simply allows us to maintain a small number of surface reactions. 

Chemical activation energy barriers, if non-zero, are also listed and referenced. Because the 
microscopic Monte Carlo simulation is very computationally demanding, we allow only H, O, CO, 
and C to accrete onto dust grains and to react, starting a chemistry that contains 29 surface 
reactions. The sticking coefficient is assumed to be unity. As the chemistry progresses and ices form 
on dust particles, the ice consists mainly of water, CO, methanol, CO2, and CH4. The formation of 
methanol is particularly interesting, because it is formed by gradual addition of H atoms to CO: CO 
-> HCO -> H 2 CO -> H3CO -»■ CH3OH. Table II shows the barriers for reactions H + CO ->• HCO 
and H + H2CO — > H3CO at different temperatures, as determined experimentally by lFuchs et al 



(J2009J). Such barriers are not as crit ical for surface reactio ns as for gas-phase processes, because 
they compete with diffusion barriers (jHerbst &: Millarll2008l ). The activation energy of the reaction 
CO + OH — > CO2 + H has not been studied on surfaces. Here we use a value, 176 K, taken from 
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Table 1: Initial Gas Phase Abundances 



Species Fractional Abundance w.r.t. kh 



Hp 


q nn y i n~ 2 

t?.UU A ±L/ 


n 2 


^ v in -1 

O A ±U 


r+ 

V 


i 9 v i n~ 4 


N 


7.6 x 10~ 5 





2.56 x 10~ 4 


s+ 


8.0 x 10~ 8 


Si+ 


8.0 x 10~ 9 


Na+ 


2.0 x 10~ 9 


Mg+ 


7.0 x 10~ 9 


Fe+ 


3.0 x 10~ 9 


P+ 


2.0 x 10- 10 


C1+ 


1.0 x 10~ 9 



Table 2: Vertical Desorption Energies 



Species E%(K) 



H 2 


430 


H 


450 


O 


800 


CO 


1150 


c 


800 


OH 


2850 


HCO 


1600 


H 2 CO 


2050 


H 3 CO 


5080 


CH 


925 


CH 2 


1050 


CH 3 


1175 


o 2 


1000 


co 2 


2575 


CH3OH 


5530 


CH 4 


1300 


H0 2 


3650 


HOOH 


5700 


H 2 


5700 
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the udfa gas-phase reaction network (http://www.udfa.net/). 



3. Monte Carlo Method 

Our Monte Carlo simulation method combines the microscopic approach for the surface chem- 



istry and the macroscopic app roach for the gas-phase chemistry (jGibson fc Bru ck 2000; IChang et al 



20051 : ICuppen Herbsti l2005). In this section, we first introduce these two methods and then show 



how they can be unified. 

To perform a microscopic Monte Carlo simulation of chemical reactions on a grain surface with 
N binding sites, we put these sites on a square lattice with dimensions LxL, where L = y/N. Events 
are executed in order of the absolute time when these events occur. There are three different types 
of events occurring on a lattice. Species accrete from the gas phase and land on random binding 
sites on the lattice. The absolute time of accretion is calculated by gas-phase kinetics, which will 
be discussed below. A species can either hop from one site to any one of its four nearest horizontal 
neighbor sites with equal probability or desorb from the surface. The hopping motion is subject to 
periodic boundary conditions on the lattice. We generate a random number X distributed uniformly 
within (0, 1) to determine whether a given species will perform hopping or thermal desorption. If 
X < b\/(b\ + 62)1 a species will hop, otherwise, it will desorb. The "waiting time" before hopping 
or desorption happens is calculated as r = — ln(X )/(b± + 62)) where X is another random number. 
This expression for r maps out a homogeneous Poisson distribution which mimics the unimolecular 
decay of a species in a given lattice point. The absolute time is then calculated as t = t + r where 
t is the time when the species first arrived at the binding site. 

Chemical reactions are results of these different types of events. For surface reactions, we 
assume that a reaction occurs instantly when two species arrive in the same site and can undergo a 
chemical reaction without activation energy. There are three mechanisms included here. The first 
one, known as the Langmuir-Hinshelwood mechanism, occurs by pure diffusion; i.e., a species hops 
to a neighboring site already occupied by another adsorbate and the two react. This mechanism is 
the major one studied by surface chemists. The second one is called the Eley-Rideal mechanism. 
Here a gas-phase species lands directly on a site where there is already a reactive adsorbate. In our 
simulations, this process is typically a minor one. We label the third one a "chain reaction" mecha- 
nism, which will be explained later. If there is activation energy, we simulate a competition between 
diffusion away from the site and tunneling under or hopping over the activation energy barrier to 



lead to reaction for both the Langmuir-Hinshelwood and Eley-Rideal mechanisms (jHerbst &: Millar 



2008). The activation energy potential is modeled as a square potential with width a = 1 A, so 
the rate of tunneling is calculated as 63 = uexp(—^\/2fj,E a ), where E a is the reaction barrier, v 
is the trial frequency, and \x is the reduced mass of reactants. Once we calculate 63, a random 
number Y uniformly distributed within (0, 1) is generated to determine if the reaction can happen 
via tunneling. If Y < 63/ (&i + 62 + 63), the reaction will happen, otherwise, the reaction will not 
happen, where we have also included the rate of thermal desorption as a competitive process. In 
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Table 3. Surface Reaction Network 



Number Reaction 



E a (K) Reference 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 



H + H - 
H + O - 
H + CO 
H + C - 
H + OH -> 
H + HCO 
H + H 2 CO 
H + H 3 CO 
H + CH -> 
H + CH 2 - 
H + CH 3 - 
O + O -> 2 
O + OH -> 2 + H 
O + HCO -»■ 
O + H3CO - 
O + CH - 
O + CH 2 
O + CH 3 
CO + OH 
C + OH - 
C + HCO -> 
CH 3 + HCO 
OH + H 2 CO 



H 2 
OH 

+ HCO 
CH 
+ H 2 
-> H 2 CO 
-»■ H3CO 
-»• CH3OH 
CH 2 
CH 3 
CH 4 



* C0 2 + H 
-> H 2 CO + OH 
CO + H 
CO + H 2 
H 2 CO + H 
C0 2 + H 
CO + H 
> CO + CH 
CH 4 + CO 
HCO + H 2 



Table 4 Fuchs et al. (2009) 



Table 4 Fuchs et al. (2009) 



OH + CH 2 -> H 2 CO + H 

C + O CO 

C + 2 -»• CO + O 

H + 2 — )■ H0 2 

H + H0 2 -»■ HOOH 

H + HOOH H 2 + OH 



176 www.udfa.net 



1200 Cuppen fc Herbst (2007) 
1400 Cuppen fc Herbst (2007) 
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addition to surface reactions, we allow H atoms to diffuse into the ice. When an H atom comes to 
a site by hopping or accretion, we determine if there is any species in the top 4 monolayers on the 
site that can react with H atoms without a chemical barrier. If there is such a species, a reaction 
without barrier happens with the topmost species within the top 4 monolayers. Otherwise, we 
check to see if there is any species in the top 2 layers on the site that can react with H atoms via 
an activation barrier. If there is such a species, we use the competition mechanism to determine if 
the reaction with a barrier can happen or not with the topmost species. 

S ome surface r eactio ns are handled in a distinctive manner. For the two surface reactions in Ta- 
ble HI Fuchs et al. (j2009l ) fit their experimental data to the hopping expression 63 = vexp(—E a /T) 
with a temperature-dependent activation energy, and a lower trial frequency of 2 x 10 11 s _1 . We 
use this frequency and these values of the activation energy in our simulation in order to match the 
experimental fitting procedure, where the lower v value is also used to calculate the competitive 
diffusion and desorption rates. 

We also must consider a phenomenon we describe loosely as a "chain" reaction. Suppose that 
there is a site with O on the topmost layer and CO underneath it. The reaction between these two 
species is not included in our small surface network because it has a sizable barrier. If a surface H 
atom diffuses onto the outermost layer of the site, it can react quickly with O to form the radical 
OH, which in turn can react with the CO molecule to form CO2 and H despite a small activation 
energy. This unusual type of "chain" reaction has already been implemented with the modified 
rate method in a three-phase gas-grain model, where the so-called "third" phase is the outermost 



layer (jGarrod &: Paulv Il201ll ). Here, because of the spatial resolution of the microscopic Monte 



Carlo simulation, we are able to treat this type of unusual reaction more rigorously. 



The macroscopic Mon t e Carlo approach is used to simulate the gas-phase chemistry (| Gibson fc Bruck 



2000 ; Vasyunin et al .112009 ; Charnley &: Rodger sl l2009h . In this approach, spatial information within 



a volume of space V surrounding a dust particle used in the simulation is not obtained. We start 
with a number of molecules for each species based on initial abundances and the volume of space. 
We then calculate the reaction rate r (s _1 ) for each reaction in the volume V. The waiting time 
for reaction i is determined by t, = — ln(Z)/rj where Z is yet another random number uniformly 



Table 4: Activation Energy Barriers 



T(K) 


CO + H 


H 2 CO + H 


12.0 


390 


415 


13.5 


435 


435 


15.0 


480 


470 


16.5 


520 


490 



Note. — The activation barriers are from lFuchs et all (|2009h 
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distributed over (0, 1). The absolute time when the i-th reaction happens is ti = tio + Ti, where tio 
is either the time when the i-th reaction happened previously or 0, when we start the simulation. 
We sort the absolute time for each reaction and find the smallest time at each step. If the i-th 
reaction happens first, at a time ti, we update the numbers of all species involved in this reaction. 
We update the rate of the k-th reaction only if the number of its reactant molecules change due to 
the i-th reaction. If the reaction rate for the k-th reaction is before the i-th reaction happens 
and r' k afterwards, the absolute time when the k-th reaction will occur again is updated to be 

t k = ti + r 4{t k -ti). (2) 
r k 

Once again, we find the reaction that will occur next, and continue the procedure. This particular 
method is, not surprisingly, known as the next reaction approach. 

The rate for the i-th gas phase reaction, where the word "reaction" includes accretion onto a 
grain, is calculated in one of three ways, as shown by the relation 

kiNj one — body reactions 

kiNkNi two — body reactions (3) 

^VT\j^£r accretion 

where Nj, Nf~, Ni are numbers of assorted reactants, ki is a rate coefficient, T g is the gas phase 
temperature, L is the dimension of the grain lattice (N 1 ^ 2 ), s is the density of sites, and m; is the 
mass of an accreting species. The next reaction method essentially models each chemical reaction 
as an inhomogeneous Poisson Process, with a rate that can change as other species change in the 
simulation. The fact that the rate is not a constant makes the merger of the two different Monte 
Carlo simulations straightforward, since desorption from the ice into the gas is just another process 
leading to a change of rate in the gas phase. 

Photon arrival events are modeled as independent homogeneous Poisson processes. We posit 
that only species in the top two monolayers can be photodesorbed, so the photodesorption rate is 
initially assumed to be 

k pd = 2F FUV Y pd -Krj, (4) 
where we assume Y p d to be 0.001 for all species, a number that is supported by experimental 



data (jOberg et al.l l2009al lbl) . The waiting time for a photodesorption is t pd = — \n{W)/k p d where 



W is a random number uniformly distributed within (0, 1). Once a photodesorption event can 
happen, we randomly assign a site where it might take place. If the site is empty, then the 
photodesorption is a null event. If there is only one layer of species on that site, we generate yet 
another random number U to determine if photodesorption events can happen or not. If U < 0.5, 
a photodesorption event occurs, otherwise, this is a null event. If there are at least two monolayers 
of species on that site, we randomly pick one of the species in the top two monolayers and allow 
this species to desorb into the gas phase. For simplicity, we ignore photodissociation in our reaction 
network. 
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The unification of the microscopic and macroscopic Monte Carlo procedures is straightforward. 
Each event in the gas-grain system has an absolute time when it happens. We sort out the smallest 
one and execute that event, which could be hopping, hopping leading to reaction, accretion, des- 
orption, or reaction in the gas phase. Accretion is considered to be a gas phase reaction. Once an 
accretion event happens, we update the number of accreting species, all necessary gas-phase reac- 
tion rates, and the absolute time of relevant reactions in the gas phase. We then add a molecule of 
the accreting species to the lattice. For a thermal desorption event, we delete the molecule from the 
lattice and add one molecule of that species to the gas phase. We once again update all necessary 
values in the gas phase. For all events other than desorption or accretion, their execution occurs 
either on the dust particle or in the gas. 

Surface proc esses in our sim ulation require more random numbers than needed with Gillespie's 



original method (lGillespielll976l ). mainly because of the extra spatial resolution capability of the 
microscopic Monte Carlo simulation. For instance, for hopping and evaporation processes, as in 
Gillespie's original method, we use one random number to determine the waiting time and then use 
the second random number to determine whether to hop or evaporate. However, in order to find the 
position of the particle after hopping, we need another random number to determine its location. 
On the other hand, it is not necessary to use the third random number for evaporation processes 
because we do not have to keep track the position of the evaporating particles. Similarly, accretion 
and photodesorption all require more random numbers to find the location when these processes 
happen. Moreover, reactions with barriers are modeled as the results of more than one encounter 
in the microscopic Monte Carlo simulation, so we need extra random numbers to determine if an 
encounter with a barrier between the two particles can lead to a reaction or not. 

Because the microscopic Monte Carlo approach is very computationally expensive, we cannot 
use a whole standard grain with 1.88 x 10 6 sites. On the other hand, previous calculat ions show that 



there is hardly any size dependence for chemical kinetics on a rough grain surface (jCuppen et al 
2009). So we only take 1% of a standard grain surface; i.e. 1.88 x 10 4 sites, consisting of a 137 x 137 



lattice. The gaseous volume associated with this part of the grain surface is V = 10 10 /n# if the 
dust-to-gas number ratio is 10 -12 . Figure [U shows that other than increased fluctuations, calculated 
abundances for selected gas phase species by a pure gas-phase Monte Carlo simulation do not vary 
with the volume of isolated gas down to a volume of V. The fluctuations shown in Figure Q] are at 
a maximum when HCO + and OH have fractional abundances of < 10~ 10 , which means that there 
is at most one molecule of each species in the volume V, so that the fluctuations are expected to be 
large. Species of higher abundance show smaller to negligible fluctuations. To improve matters for 
species of low abundance, we can always average gas abundances, as discussed in Section 4. So gas 
phase chemistry can be simulated within 1% of the volume associated with 1% of the grain without 
introducing any major error. Because 1% of a grain surface area is large enough to simulate the 
surface chemistry accurately, we can simulate the gas- grain chemistry using 1% of a grain surface 
area and an associated volume of gaseous species. 

We implemented our method using the C language and ran simulations on a Dell Precision 
T7500n Work station. The CPU running-time varies significantly depending on physical parameters 



11 



chosen. For Model 1 (see Table [5]) at 10 K, it takes only a few minutes to finish a simulation through 
2 x 10 5 yr while for Model 3, it takes about 4 hours to run an analogous simulation. Moreover, most 
of the CPU running-time is spent on late times of evolution in the system. This occurs because at 
these times H atoms have to hop many times to find partners with which to react. 



4. Results 

In this section we report the molecular abundances in the ices as they evolve on dust grains 
at different temperatures, gas-phase number densities, and ratios of diffusion barrier to desorption 
energy, p. To account for uncertain parameters, we used four different models, shown in Table [5j 
and defined by the value of the hydrogen nuclear density nn, which is either 2 x 10 4 cm 3 s _1 
or 1 x 10 5 cm 3 s _1 , and the value of p, which is either 0.77 or 0.5. The temperatures at which 
the models are run were chosen to be 10 K, 13.5 K, 15 K, and 16.5 K, because the experimental 
rea ction barriers for H + CO and H + H2CO are only available at 12 K, 13.5 K, 15 K and 16.5 



K (jFuchs et al.ll2009l ). Because 10 K is a standardly used temperature, we extended the reaction 
barriers to 10 K by assuming that the reaction rates at 10 K are the same as at 12 K as in 
Cuppen et al. (j2009l ). The dust-to-gas number ratio is fixed at 10~ 12 . We ran each simulation 
for a time period of 2 x 10 5 yr, the so-called early time for cold cores, where best agreement 
for gas-phase species is typically obtained. Moreover, the a mount of water ice produced at this 



time is sufficient to explain observa tions, as discussed below (jGibb et al.ll2004l ; iPontoppidan et al 



2004 



Boogert fe Ehrenfreundll2004l ). For each set of parameters, we ran the simulation four times 
with different random seeds and then calculated the average abundances in order to reduce the 
fluctuation in the results. 



4.1. Ice abundances at different temperatures and hydrogen densities 

We first discuss results for Models 1 and 2, where p is fixed at 0.77. For Model 1, the density of 
hydrogen nuclei is 2 x 10 cm" " 3 while for Model 2, it is 10 5 cm . Below we focus on the molecules 
in our models with major ice abundances over some temperature range. We start, however, with a 
look at the gas-phase species that accrete onto grains. 



Table 5: Parameters for Different Models 



Parameter Model 1 Model 2 Model 3 Model 4 
n H (cm" 3 ) 2.0(4) 1.0(5) 2.0(4) 1.0(5) 
E h o/E m 0.77 0.77 0.5 0.5 



Note. — a(b) means axlO 6 
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Time (yrs) 

Fig. 1. — Selected gas-phase abundances using a purely gas-phase Monte Carlo simulation with gas 
phase volumes 100 V and V at 10 K. The density of H nuclei is fixed at 2 x 10 4 cm~ 3 . 
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Figure [2] shows the gas-phase fractional abundance vs time of the four accreting species in our 
models at two different temperatures (10 K and 15 K) and nn values. Over all, we can see that the 
change of temperature does not cause much difference in the abundances of the gas phase species 
while the change in hydrogen density does. Specifically, increasing nn to 10 5 cm -3 reduces the 
abundances of all accreting species due to the larger accretion rate, especially at the latest times. 
Not all species decrease in abundance with increasing time; for example, gas-phase H and CO 
increase in abundance over the entire range of time for «h = 2xl0 4 cm -3 , while C first increases 
and then decreases at both densities, reflecting the well-known conversion of C + into CO through 
neutral atomic carbon. 

Figure E] shows the build-up with time of CO, H 2 CO, C0 2 , CH 3 OH, CH 4 , HOOH, and H 2 
in the ice using Model 1 at different temperatures. The molecular abundances on the grain mantle 
are depicted in units of monolayers, so that if we have 1 monolayer of CO in the ice, there are 
1.88 x 10 4 CO molecules, although they need not be in the same monolayer. Because the surface 
we used in the simulation is 1% of an actual grain surface for standard-sized grains, there are 
in reality 1.88 x 10 6 CO molecules present. Our dust-to-gas number ratio is fixed at 10 -12 , so 
that the total number H atoms associated with 1 grain is 10 12 ; thus, the fractional abundance of 
CO molecules on a grain surface is 1.88 x 10~ 6 . Another way of expressing the abundance of ice 
molecules is to relate them as percentages of the abundance of water, which nears 50 monolayers 
toward the end of our calculation. 

As can be seen in Figure [31 almost all of the major ice species other than hydrogen peroxide at 
early times increase in mantle abundance with time, some more sharply than others. If we look first 
at CO, we see that with increasing temperature, the CO abundance increases more strongly while 
the CH3OH abundance behaves in the opposite manner. By 16.5 K, CO is a major component, 
occupying almost 10 monolayers at the end of the simulation, while there are only 1-2 monolayers 
of CH3OH. At 10 K, on the other hand, the abundances are reversed. This result is not surprising 
because an increase in the surface temperature increases the rate of desorption of hydrogen atoms 
from the grain surface, thus reducing the H addition rate to CO. This effect is more important 
than an increase in rate due to hopping over the H + CO activation energy barrier. From Table SI 
we can see that because the activation energy rises from 390 K to 520 K as the temperature goes 
from 12.0 K to 16.5 K, the Boltzmann factor rises only slightly from 7.7 x 10 -15 to 2.1 x 10 -14 . 
The analogous change in the Boltzmann factor for the desorption rate of H is much greater since 
the desorption energy is constant. A similar statement can be made about the relative lack of 
temperature dependence for the hopping over the barrier for H + H 2 CO. Note that the abundance 
of formaldehyde shows a weak temperature dependence because both H + HCO and H + H 2 CO 
are reduced in rate due to the increasing desorption of H with increasing temperature. Our model 
can produce a moderate amount of H 2 CO between 10 K and 16.5 K. 

Model 1 produces a significant amount of CH4 at all temperatures. The desorption energies 
of C, CH, CH 2 , and CH3 are sufficiently high that these species will not desorb rapidly through 
16.5 K, but will be converted to CH4 with high efficiency. The major alternative is reaction with O 
atoms to form CO. These processes are not the dominant ones for CO except at very early times, 
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Fig. 2. — The evolution of the abundances of gas-phase species that accrete onto grain surfaces as 
followed with Models 1 and 2. Panel (a): T = 10 K, Model 1, panel (b): T = 10 K, Model 2, panel 
(c): T = 15 K, Model 1, panel (d): T = 15 K, Model 2. 
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Fig. 3. — The abundances of water, hydrogen peroxide, and major carbon-containing species in the 
ice mantle are shown in terms of monolayers as a function of time at different temperatures using 
Model 1. 
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Fig. 4. — The abundances of water, hydrogen peroxide and major carbon-containing species in the 
ice mantle are shown in terms of monolayers as a function of time at different temperatures using 
Model 2. 
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after which accretion of CO from the gas dominates. Carbon dioxide is a major carbon-bearing 
molecule that can be formed between two heavy species on grains between 10 K and 16.5 K. This 
molecule is produced mainly by the reaction OH + CO — > CO2 + H, which occurs by an unusual 
"chain mechanism" discussed earlier, in which an H atom first combines with an O atom lying above 
a CO molecule, so that the OH need not undergo horizontal diffusion to react with CO. Water, 
the dominant species, is formed by a variety of surface reaction sequences starting with atomic O, 
molecular O2, and even O3. Use of only the first two syntheses (see Table [3j) is sufficient in the 
temperature range studied to pr oduce 30 monolayers in 2 x 10 5 yr, far more than could be produced 



by ac cretion of gas-phase water (jGibb et alj|2004l ; iPontoppidan et al.ll2004l ; iBoogert &: Ehrenfreund 



2004). The O + H mechanism is dominant at 10 K while O2 + H is the dominant mechanism at 
16.5 K. There are also moderate amounts of HOOH produced on grain surfaces in Model 1. The 
production of hydrogen peroxide generally decreases with higher temperature due to an increase in 
rate for the reaction H + HOOH — > OH + H2O, which possesses a barrier. We can also see that 
the abundance of HOOH has a bump at 16.5 K at early times. This can be explained by our initial 
abundances and the rapid increase of H atoms in the gas phase at early times. 

Analogous results with Model 2 are shown in Figure dj which indicates the influence of higher 
density. One obvious consequence of the increased density is the increased total thickness of ice 
on the grain surface, which occurs because of the increased deposition flux. For example, at 
temperatures above 10 K, 50 monolayers of water ice can be produced in 2 x 10 5 yr. The other 
major difference is that greater amounts of CO and H2CO are produced on grain surfaces, especially 
at 10 K. This increase in the mantle CO abundance occurs for two reasons. First, the fractional 
abundance of gas-phase atomic H in the higher density case decreases strongly after 2 x 10 yr, as 
shown in Figure such that it is far lower than that of CO. The lowered relative flux of H onto 
grains compared with CO makes it more difficult to destroy CO. Second, CO molecules are more 
easily buried before H atoms can react with them with the increased flux of species. 

The increased abundance of H2CO is more difficult to understand because it is more difficult 
to convert CO to H2CO; however, it is also more difficult to convert H2CO to CH3OH at higher 
density. Based on our simulation, the absolute abundance of H2CO increases, but the abundance 
ratio of H2CO to CO does not increase at higher density. The increased abundance of HOOH can 
be understood in similar way; i.e., the rate of the formation reaction (H + HO2) decreases less 
rapidly than the rate of the destruction reaction (H + HOOH). From Figure HI we can see that the 
absolute abundance of CH3OH is somewhat decreased at higher density, while the abundance ratio 
of CH3OH to CO clearly drops precipitously, as explained earlier. 

Table [6] and Table [7J show the abundances of water ice, hydrogen peroxide, and all major 
carbon-bearing species in the mantle ice at the final time of 2 x 10 5 yr. Table [6] lists results for 
Model 1, while Table [7J lists results for Model 2. The abundance of water is given as its fractional 
abundance with respect to the gas-phase nuclear hydrogen density nn, while the abundances of the 
other ice species are given as the percentages relative to water ice. Some large differences can be 
seen between results at the two different densities, especially at low temperatures. For example, 
the abundance of CO with respect to water is a factor of 6 lower at low density at 10 K, and only 
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Table 6: Model 1 Abundances of Major Mantle Species 



Species 


10 K 


13.5 K 


15 K 


16.5 K 


CO 


5.7 


8.4 


12.8 


26.8 


H 2 CO 


4.6 


4.0 


4.2 


3.6 


co 2 


19.6 


13.4 


15.2 


21.4 


CH3OH 


27.7 


17.3 


10.7 


5.4 


CH 4 


18.5 


20.6 


20.9 


15.0 


HOOH 


4.8 


0.77 


0.44 


0.66 


H 2 


4.7(-5) 


6.4(-5) 


6.6(-5) 


6.0(-5) 



Note. — The abundances pertain to a time of 2 x 10 5 yr. 

Note. — The abundance of water ice is the fractional abundance with respect to gas- phase i%h, while the abundances 
of the other ice components are percentages with respect to the water ice. a(-b) means axlCT 6 . 



Table 7: Model 2 Abundances of Major Mantle Species 



Species 


10 K 


13.5 K 


15 K 


16.5 K 


CO 


34.9 


48.1 


53.5 


58.5 


H 2 CO 


22.4 


9.7 


6.8 


4.4 


co 2 


77.3 


38.7 


29.5 


26.6 


CH30H 


27.2 


8.9 


5.0 


2.4 


CH 4 


15.5 


12.0 


10.7 


8.4 


HOOH 


19.0 


7.1 


4.8 


2.8 


H 2 


5.6(-5) 


8.9(-5) 


9.9(-5) 


1.0(-4) 



Note. — The abundances pertain to a time of 2 x 10 5 yr. 

Note. — The abundance of water ice is the fractional abundance with respect to gas- phase rin, while the abundances 
of the other ice components are percentages with respect to the water ice. a(-b) means axlCT 6 . 
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a factor of 2.2 lower at 16.5 K. These results will be compared with observations in Section [5j 

Now let us consider the distribution of assorted species throughout the grain mantle at 2 x 10 5 
yr. Note that the inner layers were formed at earlier times and remain relatively passive. Figure [5] 
shows the fraction of each monolayer occupied by carbon-bearing species, hydrogen peroxide, and 
water starting on the surface of the grain and ending at the uppermost and partially occupied layers 
for Models 1 and 2 at 10 K and 15 K. Because the gas-phase atomic C fractional abundance peaks 
at early stages, its accretion onto dust particles and conversion to methane by successive reactions 
with atomic hydrogen means that the abundance of CH4 should peak in the inner layers of the ice. 
Indeed, in Figure EJ CH4 generally occupies follows this predicted trend, except for a secondary 
peak at higher monolayers at 15 K for the higher density. This anomaly can be explained by the 
depletion of O, which is a competitor of atomic hydrogen in reacting with C and carbon hydrides. 

The abundance of CO is typically relatively flat as a function of monolayer until the outermost 
layers, but actually peaks strongly in the innermost layers for Model 1 at 10 K. This result is 
non-intuitive because CO n ormally forms on g rains by accretion from the gas, and its production 



in the gas takes some time (lOberg et al.ll201ll ). There are two explanations for a relatively large 
amount of CO at inner monolayers in our simulation. First, CO can be produced on grain surfaces 
by the reaction C + O — > CO. There are high abundances of gas phase O and C at very early times 
of our simulation; thus, upon accretion of these atoms onto the ice, a significant amount of CO can 
be produced. Secondly, the abundance of gas phase H atoms is small at early stages, so the flux 
of H atoms accreting onto the grain surface is small, which reduces a standard destruction route 
for CO ice. To confirm this hypothesis, we increased the initial fractional abundance of atomic 
H to 4 x 10 -4 , and found that the CO abundance at the inner layers virtually disappears. At 
higher densities, more CO is produced at later times by accretion, while at 15 K and lower density, 
it is harder to destroy CO with H atoms due to the increase in desorption. The CO abundance 
distribution is also discussed in Section 5, where we compare our results with observations and 
inferences. 

The methanol abundance typically peaks at late stages, especially when the gas phase CO is 
large. The rapid increase of methanol at outer layers at 15 K for higher densities is particularly 
interesting. This effect can be explained by the heavy depletion of gas phase O under these condi- 
tions, which leads to a lower accretion flux, so that H atoms on the grain are more likely to convert 
CO into methanol than to react with surface O. 

The abundances of CO2 and water at different layers are close to constant except towards 
the outermost layers. The abundance of H2CO at different layers behaves roughly in the same 
manner as CO regardless of nn- The distribution of hydrogen peroxide is dependent on the gas 
phase density. At lower densities, HOOH drops quickly from inner to outer layers. The effect at 
the outer layers is explained by the rapid increase of the gas phase ratio of H to O, and the fact 
that HOOH can be converted to water efficiently in our temperature range. At higher densities, 
HOOH maintains a moderate abundance at the inner layers before dropping rapidly at the outer 
layers. The effect at the inner layers can be understood by a drop of gas phase H and the fact that 
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HOOH can be buried more easily at higher densities. 

Finally, we show a picture of the icy surface for Model 1 at 2 x 10 5 yr at 10 K in Figure El The 
figure shows the number of layers at each site on the grain surface. We can see that the surface 
has actually grown to be very rough, so that our explanations should be taken with the proverbial 
grain of salt. 

4.2. Influence of the ratio of diffusion barrier to desorption energy 

In this subsection, we discuss results when p is set to a lowered value of 0.5, which increases 
the diffusion rate for a fixed desorption energy. We define Model 3 and Model 4 in analogy with 
Models 1 and 2 except for the lower value of p (see Tabled]). It is tedious to finish a microscopic 
Monte Carlo simulation with the faster diffusion rate at temperatures higher than 10 K, so we only 
give results at 10 K. 

The temporal evolution of accreting gas phase species is very much like that in Figure (2j as 
expected. Figure [7] shows the temporal evolution of water, hydrogen peroxide, and major carbon- 
bearing ice species on dust particles with Models 3 and 4, while Figure [8] shows the fraction of each 
monolayer occupied by these species at the final time of 2 x 10 5 yr for these models. Compared 
with results from the slower diffusion rate (p = 0.77) at the same temperature, one major difference 
is that more CO and less CH3OH are produced. The increase in the rate of diffusion affects the 
competition between diffusion and reaction for Langmuir-Hinshelwood-type reactions with chemical 
activation energy barriers, for example H + CO and H + H2CO. If the diffusion rate increases, it is 
more likely that H atoms will diffuse away from a site than react. On the other hand, faster diffusion 
can increase the probability that an H arrives to encounter a CO molecule at a site. Under the 
conditions of the simulation, the former effect appears to dominate. Table [8] shows the abundances 
of all major species on a dust grain at 2 x 10 5 yr compared with the water abundance, as in Tables [6] 
and[7l We see that with the faster diffusion rate, increasing «h can strongly increase the percentage 
of CO molecules on the grain surface as also occurs with the slower diffusion rate. If we consider 
the total absolute abundance at 10 K of the major carbon-bearing species, we see that raising the 
diffusion rate does not change this value because the decrease in total abundance with respect to 
water ice is balanced by an increase in the water ice abundance (see Table [8] as compared with 
Tables [6] and [7|) . With a faster diffusion rate, the temporal evolution of HOOH and the distribution 
of HOOH at different layers behave similarly to the p = 0.77 results at higher temperatures, which 
is reasonable because the major effect of increasing the temperature is to increase the diffusion rate 
of species. Such an analysis is useful for other species as well. 

5. Comparison with Observation and Other Models 

The primary purpose of this paper is to introduce a stochastic method to merge two different 
types of Monte Carlo simulations. In doing so, we were forced to confine ourselves to a rather small 
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surface reaction network with only 29 reactions, which is mainly useful for calculating abundances 
and environments of water ice, hydrogen peroxide, and the major carbon-containing ices. Never- 
theless, it is still interesting to compare our calculated abundances for major ice species with recent 
observations and results from other surface models. Comparison with gas-phase species shows that 
for the time interval calculated, very little difference exists between the results here and results of 
gas-phase models. For comparison with observed ice abundances, we utilize our results at a time 
of 2 x 10 5 yr. 
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Fig. 5. — The fraction of each monolayer occupied by water, hydrogen peroxide, and major carbon- 
bearing species using Models 1 and 2 shown as a function of monolayer at the final time of 2 x 10 5 yr. 
Panel (a): T = 10 K, Model 1, panel (b): T = 10 K, Model 2, panel (c): T = 15 K, Model 1, panel 
(d): T = 15 K, Model 2. 
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Grain Base 



Fig. 6. — Number of layers at each site on the grain surface after 2 x 10 5 yr of evolution at 10 K 
for Model 1. 
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Fig. 7. — The abundances of assorted species in the ice mantle are shown in terms of monolayers 
as a function of time at 10 K using Model 3 (Panel (a)) and Model 4 (Panel (b)). 



25 



(a) 



(b) 









10 



20 



30 



40 



50 



60 
Layer 



20 



40 



60 



80 



100 



Fig. 8. — The fraction of each monolayer occupied by water, hydrogen peroxide, and major 
carbon-bearing species using Models 3 (Panel (a)) and 4 (Panel (b)) at 10 K shown as a function 
of monolayer at the final time of 2 x 10 5 yr. 



Table 8: Abundances of Major Mantle Species at 10 K 
Species Model 3 Model 4 



CO 


10.5 


48.7 


H 2 CO 


8.5 


10.9 


co 2 


11.3 


28.6 


CH 3 OH 


9.8 


4.1 


CH 4 


21.2 


10.1 


HOOH 


0.50 


4.5 


H 2 


5.9(-5) 


9.6(-5) 



Note. — The abundances pertain to a time of 2 x 10 5 yr. 

Note. — The abundance of water ice is the fractional abundance with respect to gas-phase nu, while the abundances 
of the other ice components are percentages with respect to the water ice. a(-b) means axlO -6 . 
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Table 9. Comparison of Observational and Theoretical Ice Abundances 



CO C0 2 CH3OH CH 4 HOOH 



Our Models 


Model 1, 10 K 


5.7 


19.6 


27.7 


18.5 


4.8 


Model 1, 13.5 K 


8.4 


13.4 


17.3 


20.6 


0.77 


Model 1, 15 K 


12.8 


15.2 


10.7 


20.9 


0.44 


Model 1, 16.5 K 


26.8 


21.4 


5.4 


15.0 


0.66 


Model 2, 10 K 


34.9 


77.3 


27.2 


15.5 


19.0 


Model 2, 13.5 K 


48.1 


38.7 


8.9 


12.0 


7.1 


Model 2, 15 K 


53.5 


29.5 


5.0 


10.7 


4.8 


Model 2, 16.5 K 


58.5 


26.6 


2.4 


8.4 


2.8 


Model 3, 10 K 


10.5 


11.3 


9.8 


21.2 


0.50 


Model 4, 10 K 


48.7 


28.6 


4.1 


10.4 


4.5 


Previous Models 


Master equation, 10 K 


2.6(-6) 


0.025 


35.7 


20.1 




Master equation, 20 K 


30 


0.05 


13 


1.4 




RE, slow diffusion 10 K 


27.9 


8.4(-14) 


8.4(-10) 


23.7 


4.1(-6) 


RE, slow diffusion, 15 K 


57.5 


1.3(-7) 


5.1(-5) 


5.6 


5.7(-3) 


RE, fast diffusion, 10 K 


36.3 


2.4(-7) 


0.63 


14.3 


1.8(-8) 


MRE, slow diffusion 10 K 


29.0 


8.6(-13) 


1.3(-9) 


24.7 


7.3(-6) 


MRE, slow diffusion, 15 K 


54.8 


1.4(-5) 


l.l(-4) 


7.2 


3.9(-3) 


MRE, fast diffusion, 10 K 


31.1 


1.9(-6) 


2.0 


14.3 


3.4(-8) 


Three phase model, 10 K 


10 


5.7 


1.1 


30 




Observations 


Toward Low Mass Stars 


29 


29 


3 


5 




Toward High Mass Stars 


13 


13 


4 


2 




Toward Background Stars 


31 


38 


4 






Elias 29 










<13 


Mon R2/IRS3 










<5 


Mon R2/IRS2 










<13 


W51/IRS2 










<16 



Note. - The abundances of the ice components are percentages with 
respect to the water ice. a(-b) means axl0~ b . Observational results for 
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carbon-bearing species are from the compilation of Obcre et al. 



servational results for hydro gen peroxide are from ISmith et al.l (120111 ). mas 



(1201 ob- 



ter equation results are fromlStantcheva &; Herbstl (|2004l ) , while three-phase 



model results are from iGarrod &: Paulv I ( 2011 
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The median fractional abundance of water ice in dense cold sources with respect to »h is about 



5 x 10" 5 (Gibbetal 



2004 



Pontoppidan et al. 



2004 



Boogert &: Ehrenfreund 



2004 



Oberg et al 



201 ll ). which is in ge neral agreem e nt wi th our results for all models. A compilation of Spitzer and 
ISO observations by lOberg et al.1 (120111 ) shows that the median abundances of CO, CO2, CH3OH 
and CH4 ice respectively as percentages of water ice are 29, 29, 3, 5 toward low mass stars, 13, 13, 
4, 2 toward high mass stars, and 31, 38, 4, - towards bac kground star s , resp ectively. In Table we 
summarize the ice abundances from the observations of lOberg et al.l (|201ll ) along with the results 
of assorted simulations in this paper and others. Our simulations differ in their agreement with the 
observed results, although it is difficult to use the comparison to constrain both the temperature 
and the density. In general, our methane results are significantly too high, while simul ation results 



for th e other three carbon species are much more reasonable. For example, values of lOberg et al 



(|201ll ) towards low-mass stars are quite reasonably fit by our Model 1 and Model 2 results at 16.5 K 
while the observations towards high-mass stars are fit reasonably by Model 3 results at 10 K, with 
the exception of methane. The temperature dependence is non-intuitive because one might expect 
that the ices in front of high-mass and luminous sources possess a higher temperature than ices in 
front of low-mass and background stars. Significantly higher ice temperatures wou ld lead to less 
of the volatile species, such as CO and CH4, in the form of ices (jOberg et al.ll201ll ). Our models, 



on the other hand, associate a low CO ice abundance with a lower temperature, at which CO can 
be converted chemically to methanol. One must remember that with the fixed temperatures used 
in our model, we cannot distinguish between the temperature at which ice species form and the 
temperature, which may differ, at which they are observed. Thus, the predicted ice abundances in 
our model are most appropriately compared with observations with background stars. The general 
degree of agreement can be seen in the histogram displayed in Figure (H where we have also included 
Model 4 at 10 K. 

We can also compare our results with respect to water ice with the results of other simulations, 
as shown in Table [9l Using an approximate master equation stochastic approach over a time of 10 



vr. iStantcheva &: Herbstl (|2004l ) calculated the abundances of CO, CO2, CH3OH and CH4 ice as a 
percentage of water ice. Their reported results at 10 K and 20 K at a time of 3.2 x 10 5 yr show very 
low abundances of CO2 ice at both temperatures, and almost no CO at all at 10 K. Instead, the 
conversion to methanol is much too efficient because they used an ultra-fast diffusion rate {E^/E-q 
= 0.3). On the positive side, the master equation results do show a low methane abundance at 20 
K. We also used our network to calculate results with both rate equation (RE) and modified rate 
equation (MRE) approaches to the surface chemistry with a so-called "two-phase" approach. In this 
approach, no distinction is made between the surface and interior monolayers. The RE and MRE 
simulations were undertaken at a density of 2 x 10 4 cm" 3 for both slow and fast diffusion. We see a 
strong similarity between the results at 2 x 10 5 yr and the results of the master equation approach, 
where there is hardly any CO2 formed. Moreover, the production of CH3OH is not efficient enough 
except for the MRE model with fast diffusion (E^/Ej) = 0.5) at 10 K. Recently, three-phase models, 
in which the surface is distinguishe d from the interior monolayers, have been un dertaken with the 



2011 



rate e quation GRAINOBLE Model (ITaauet et al.ll2012r). modified rate equation s (jGarrod h Pauly 



) and macroscopic Monte Carlo approaches (|Charnley fc Rodgerd 12009). Here, we list the 
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results for the simulation by iGarrod fe Paulv I (|201ll ). which partially includes stochastic effects 



at 10 K and 2 x 10 5 yr. This simulatio n is more success f ul exc ept for the high abundance of 



methane. Both the three-phase model of IGarrod Pauly I (|201ll ) and our simulations achieve a 
reasonable abundance of CO2 at least partially because of the "chain reaction" mechanism discussed 
earlier. Our CO2 c alculated abundances are in somewhat better agreement with the observation of 
Oberg et al.l fcoilh . 



Recent analysis of spectrosco pic data from diffe rent sources gives upper limits to the HOOH 

L). We can see in Table [9] that most of our 



abundance relative to water ice (Smith et al 
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simulation results are below these upper limits, although there is far more HOOH produced in our 
models than by rate equation or modified rate equation approaches. On the other hand, recent 
observational studies show that the fractional abundance of gas phase HOOH toward p Oph is 



about 10 (iBergman et al. 1120111 ) while in Model 1 at 10 K, the gas phase HOOH fractional 



abundance is about 4 x 10~ 10 . This shows that moderate coverage of HOOH, about 1% of the 
water abundance, is necessary to explain the gas phase abundance of HOOH unless a more efficient 
desorption mechanism exists in the system. 

It is also interesting to compare the distribution among the monolayers of major carbon-bearing 
specie s obtained from ou r simulation at 2 x 10 5 yr with inferences from observa tions and laborator y 
work ( Oberg et al. 2011 ). by which one ca n distinguish amon g environments ( Tielens et al.l[l99ll ). 
Although the situation is rather complex, lOberg et al.l (|201ll ) state that "In general, observed ice 
features can, however, be explained by an early ice chemistry phase characteristic of hydrogenation 
of atoms, followed by atom-addition reactions in a CO-dominated ice phase." Moreover, it appears 
that the monolayer distributions of CO and CO2 are anti-correlated; the fraction of CO found in 
a CO2 environment and the fraction of CO2 found in a CO environment are not high. The normal 
interpretation is that the CO formed early via accretion and surface reactions is mainly converted 
into CO2 (probably via CO + OH). This conversion becomes inefficient at later times, where, for 
the low mass YSO case, there can be a high abundance of almost pure CO, with some conversion 
to formaldehyde and methanol. To compare these inferences with our theoretical results, we focus 
on the final time of our calculations (2 x 10 5 yr) at which time the water ice abundance has become 
large enough to agree with observations. 

Rather than trying to use our detailed figures, we divide the number of monolayers of ice into 
three equal groups - inner, middle, and outer, and list in Table [10] the fraction of the monolayers 
in each group occupied by assorted molecular ices. Of the different model results in this table, 
we see that Model 2 at 15 K appears to agree with the observational data best. In particular, 
methane decreases from a fractional occupancy of 0.096 in the inner monolayers to an occupancy 
of 0.024 in the outer monolayers, while CO2 decreases more weakly from 0.15 (inner monolayers) 
to 0.13 (outer monolayers). On the other hand, both CO and its hydrogenation product, CH3OH, 
increase in abundance from inner to outer monolayers. Thus, we conclude that methane and carbon 
dioxide tend to be produced preferentially at early times, that CO and CH3OH tend to be produced 
preferentially at later times, and that CO and CO2 are in a sense anti-correlated spatially. Other 
models fit the observational inferences more poorly, although Model 4 at 10 K is also reasonably 
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good. If we wish to divide environments into those that are polar and those that are non-polar, we 
can compare the sum of the abundances of non-polar (or relatively non-polar ) species - CO, CO2, 
O2 — with the sum of the abundances of polar species - H2O, HOOH, H2CO, and CH3OH. If we do 
this, we find that there is no environment in which the non-polar species are in the strong majority. 
If there is a decidedly non-polar environment, it would have to occur within a smaller group of 
monolayers than those groups considered here. Such an environment appears to occur at very high 
monolayers in Model 2 at 15 K and Model 4 at 10 K, as shown in Figures [5] and [81 where it can 
be seen that the CO and water abundances cross near the uppermost monolayers. Nevertheless, it 
is fair to state that the model does not really reproduce the existence of a non-polar environment 
where CO freeze-out is dominant. 
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Table 10. Distribution of Abundances into Inner (I), Middle (M), and Outer (O) Environments 
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Fig. 9. — Histograms of the abundances of major carbon-bearing surfa ce species by s imula tion and 



observation with respect to water ice. Observational results are from lOberg et al.l (|201! 
1 and Model 2 results are at 16.5 K while Model 3 and Model 4 results are at 10 K. 
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6. Discussion and Conclusions 



For the first time, we have performed a unified microscopic-macroscopic Monte Carlo simulation 
of gas-grain chemical modeling of a cold dense cloud. In this approach, the gas phase chemistry, 
which consists of more than 4000 reactions, is modeled by a macroscopic Monte Carlo simulation, 
while the surface/mantle chemistry, which consists of 29 reactions, is modeled by a microscopic 
Monte Carlo simulation, in which the position of each individual species is followed. For simplicity, 
we only allow H, O, C, and CO to accrete onto grain surfaces. Although our surface reaction 
network is strongly simplified compared with a full surface reaction network, we can produce both 
water and the major carbon-bearing species on dust grains for a cold interstellar core. 

The unified simulation is made possible by modeling gas phase reactions as inhomogeneous 
Poisson processes and adopting the "next reaction" simulation method. In this approach, each 
chemical reaction in the gas phase has an absolute time for the reaction to start, which will change 
due to the occurrences of other chemical reactions in the gas phase that change the abundances of 
its reactants. The advantage of modeling gas phase reactions in this manner is that any desorption 
processes from grain surfaces can be treated in the same way as gas phase chemical reactions, thus 
making it facile to couple gas phase proc esses with grain surface processes. This approach is superior 
to Gillespie's algorithm (|Gillespidll976l ). which models gas phase reactions as homogeneous Poisson 
processes, in which a bundances are not changed by other processes. Our current model is superior 
to our earlier model (jChang et al.l 120071 ). which couples a microscopic Monte Carlo simulation for 
the grain surface with rate equations for the gas, because the restrictions required in our earlier 
model have been removed. In addition, most of the CPU time used in both approaches is spent on 
the microscopic Monte Carlo simulation of surface chemistry, so our new simulation does not lose 
much efficiency although the gas-phase chemistry is simulated by a Monte Carlo method instead 
of the more efficient rate equation approach. 

Our results for gas-phase abundances show few differences from results obtained from rate 
equations, which means that they are in good agreement with observation at the time of 2 x 10 5 yr 
used for the surface chemistry. This result is not surprising; IStantcheva &: Herbstl ([2004]) showed 
that differences can show up only at longer times, where results on surfaces begin to affect gas-phase 
results more strongly. Although our network does not contain surface reacti ons to produce comple x 
organic molecules in cold cores, as recently detected in the gas phase by iBacmann et al.l (|2012l ). 
we do produce gaseous methanol (CH3OH) by photodesorption from the ice. With Model 1 at 10 
K, we obtain a fractional abundance for methanol in the gas of ~ 6 x 10~ 10 at the final time of 
the c alculation, which agrees well with the observed abundance in TMC-1 of ~ 10~ 9 ( Smith et al 
20041 ). 



Our results for the icy mantles can be compared with observations of total mantle abundances 
of assorted ices as well as their local environments in the ice. Regarding total mantle abundances, 
our results show that after a time of 2 x 10 5 yr, a number of models with different temperatures, 
densities and diffusion barrier-to-desorption energies can reasonably reproduce the observed abun- 
dance of water ice as well as the main carbon-bearing species in front of high-mass YSO's, low-mass 
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YSO's, and backgr ound stars, with the exception of methane, which is significantly overproduced 



(jQberg et al.l 120111 ). Regarding local environments, if we divide the total number of monolayers 
into three groups (inner, middle, and outer) we can reproduce aspects of the inferred monolayer 
distributions and correlations for various species, but there is no one model that fits all of the 
information inferred from observations with high mass, low mass, and background star sources. 

Although our model can successfully couple the time-dependent evolution of gas phase chem- 
istry and grain surface chemistry, it is still not complete. First, the surface reaction network is 
still too small. A complete surface reaction network typically consists of at least a few hundred 
reactions; however, it is not clear if our model can treat this more complicated case without more 
advanced computer techniques and/or approximations. Furthermore, to model sources where there 
is a change in the physical conditions as the chemistry occurs, such as a collapse or a warm-up, 
new algorithms will have to be developed to couple these changes into the unified Monte Carlo 
approach. 
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